Event-by-event gluon multiplicity, energy density and eccentricities at RHIC and LHC 
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The event-by-event multiplicity distribution, the energy densities and energy density weighted 
eccentricity moments e n (up to n = 6) at early times in heavy-ion collisions at both RHIC 
(^/s = 200 GeV) and LHC (y^s = 2.76 TeV) are computed in the IP-Glasma model. This framework 
combines the impact parameter dependent saturation model (IP-Sat) for nucleon parton distribu- 
tions (constrained by HERA deeply inelastic scattering data) with an event-by-event classical Yang- 
Mills description of early-time gluon fields in heavy-ion collisions. The model produces multiplicity 
distributions that are convolutions of negative binomial distributions without further assumptions or 
parameters. The eccentricity moments are compared to the MC-KLN model; a noteworthy feature 
is that fluctuation dominated odd moments are consistently larger than in the MC-KLN model. 



I. INTRODUCTION 



Event-by-event relativistic hydrodynamic models do a 
good job of describing the bulk features of heavy-ion 
collisions at the Relativistic Heavy-Ion Collider (RHIC) 
and the Large Hadron Collider (LHC) However, the 
hydrodynamic description of the collective properties of 
this matter is sensitive to details of the initial conditions, 
the spatial distributions of "frozen gluon states" within 
the incoming nucleons, the event- by- event fluctuations 
in the number and energy distributions of produced par- 
tons, and not least, the thermalization of the initial non- 
equilibrium Glasma matter formed shortly after the col- 
lision. The significant amount of data now available on 
moments of flow distributions and particle spectra offer 
the promise of being able to learn about these interesting 
aspects of collective Quantum Chromo-Dynamics (QCD) 
as well as the possibility of disentangling their properties 
from the hydrodynamic transport properties of a ther- 
malizcd Quark-Gluon Plasma. 



In a previous letter Q, we introduced the IP-Glasma 
model of early time dynamics which combines the IP- 
Sat (Impact Parameter Saturation Model) model 0, H[ 
of high energy nucleon (and nuclear) wavefunctions with 
the classical Yang-Mills (CYM) dynamics of the Glasma 
fields produced after the heavy-ion collision [5-9]. The 
IP-Sat model is formally similar to the classical Color 
Glass Condensate (CGC) McLerran- Venugopalan (MV) 
model for dipole cross-sections for nucleons and nu- 
clei [To| but additionally includes Bjorken x and impact 
parameter dependence through cikonalized gluon distri- 
butions of the proton that are constrained by HERA 
inclusive and diffractive e+p deeply inelastic scatter- 
ing (DIS) data ll| and by the available nuclear fixed 
target data from EMC and E665 experiments [f2[. A 
key dimensionful quantity which determines the quali- 
tative behavior of cross-sections is the saturation scale 
Q s (x,b) which is determined sclf-consistently from the 
dipole cross-sections. The model yields good %-squarcd 



fits3 to available c+p small x HERA collider data [TT| 
and to fixed target e+ A DIS data [l2| ■ 

The IP-Glasma model relates the DIS constrained nu- 
clear dipole cross-sections to the initial classical dynamics 
of high occupancy gluon "Glasma" fields after the nuclear 
collision. Given an initial distribution of color charges in 
the high energy nuclear wavefunctions, the IP-Glasma 
framework computes the strong early time multiple scat- 
tering of gluon fields by event-by-event solutions of Yang- 
Mills equations. It is more general and accurate relative 
to ll k± factorization" models often used in the small x 
literature, especially for the soft dynamics that charac- 
terizes hydrodynamic flow. 

The IP-Glasma model has a significant feature that is 
of extreme importance for event-by-event flow studies. It 
naturally includes the effect of several sources of quan- 
tum fluctuations that can influence hydrodynamic flo\\0. 
An important source of fluctuations, generic to all models 
of quantum fluctuations, are fluctuating distributions of 
nucleons in the nuclear wavefunctions. In addition there 
are fluctuations in the color charge distributions inside 
a nucleon. This, combined with Lorcntz contraction, re- 
sults in "lumpy" transverse projections of color charge 
configurations that vary event to event. The scale of 
this lumpiness is given on average by the nuclear satura- 
tion scale Q s which corresponds to distance scales smaller 
than the nucleon size [12J. For each such configuration 
of color charges, the QCD parton model predicts dynam- 
ical event- by- event fluctuations in the multiplicities, the 
impact parameters and the rapidities of produced glu- 
ons [Tl. 



1 Since the work of Ref. [T3 |. further data and combined analy- 
ses of the ZEUS and HI collaboration results is available. It is 
important to revisit the IP-Sat fits to take this information into 
account. Work in this direction is in progress but is outside the 
scope of the present work. For nuclei, data from d + Au colli- 
sions as well as future p+A and e+A collider experiments will 
help significantly constrain Q B . 

2 A qualitatively similar model of the effect of initial state Glasma 
fluctuations (on the scale 1/Q S ) on hydrodynamic flow moments 
can be found in Ref. [T3 |. 
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These dynamical fluctuations (to be distinguished from 
geometrical fluctuations in the nucleon positions) are 
captured in the CGC Glasma flux tube picture of the 
heavy-ion collision. In this picture, the gauge fields pro- 
duced in each event after the nuclear collision are ob- 
tained by solving the CYM equations for given color 
charge densities in the two incoming nuclei and are there- 
fore functionals of these densities. The n-particlc gluon 
probability distribution can be extracted from n-particle 
gluon cumulants obtained by taking appropriate prod- 
ucts of the gauge fields and averaging these over the color 
charge density weight functionals that are given by the 
CGC effective theory. An equivalent method (followed 
here) for generating the n-gluon multiplicity distribu- 
tion is to compute single inclusive multiplicities event- 
by- event. 

In the perturbative regime Q s <C fcj_, it was shown 
that the n-particle multiplicity distribution is described 
by a negative binomial distribution characterized by the 
mean multiplicity and a parameter k [l5j]. This param- 
eter, which is unity for a Bosc-Einstcin distribution and 
infinity for a Poisson distribution, is computed to be 

k = C~2^~QlS±, where S \ is the transverse overlap area 
of the collision. In Ref. [la ], it was shown from a non- 
perturbative computation of the double inclusive gluon 
distribution in the Glasma that the Glasma flux tube 
picture is robust. Fits based on k± factorization to P+P 
and A+ A multiplicity distributions give ( ~ 1 /6 fl7l.ll8j]. 
A value of C of this order also appears to be required 
for a reasonable estimate of the amplitude of the long 
range di-hadron "ridge" correlation in p+p collisions [l!| . 
The event-by-event solutions of CYM equations, as will 
be discussed in this paper, will allow us to first obtain 
the important result that the NBD distribution is non- 
perturbatively robust, and to determine the parameters 
k and the mean multiplicities n. 

With the initial conditions generated from the IP-Sat 
model, the solution of the classical Yang-Mills (CYM) 
equations in each event allows one to determine the evo- 
lution of the spatial distributions of the produced matter 
as a function of proper time. These are characterized by 
the spatial eccentricity moments which can be defined as 

a/ (r n cos(n(f>)) 2 + (r n sin(n0)) 2 

where (•) is the energy density weighted average. These 
spatial eccentricity moments are in turn converted to mo- 
mentum space anisotropics by hydrodynamic flow. How 
efficiently this is done can in detail be determined from 
the harmonic flow coefficients v n defined through the ex- 
pansion of the azimuthal particle distribution as 

f = ^(l + E(^cos(n^. (2) 

Thus if viscous hydrodynamics is applicable, and the 
equation of state of the hot Quark-Gluon Plasma un- 
der control, one can in principle extract information on 



the transport coefficients of the flow by determining ab 
initio the initial spatial eccentricity moments e„ on the 
one hand and measuring the anisotropy moments v n on 
the other. Varying the centrality of the collision, the nu- 
clear species available, and the energy of the collision, 
likely provides sufficient handles to extract the shear and 
bulk viscosities of hot QCD matter with some degree of 
confidence. 

The IP-Glasma model can be compared and contrasted 
with previously conceived models of the initial conditions 
which include various Monte Carlo based realizations of 
the Glauber model (see (20j). and CGC inspired models 
such as KLN [UHj], f(actorized)KLN [23[, and rcBK 
[24| models. The IP-Glasma model follows all these mod- 
els by sampling nucleon positions stochastically from a 
Woods-Saxon distribution. The nucleons' color charge 
distributions are constrained on the average to reproduce 
HERA data that provide information on the spatial and 
energy dependence of gluon distributions in the proton. 
These color charge densitities are added at every trans- 
verse position and the total nuclear color charge den- 
sity is sampled to produce the color charge distribution 
in a single event. Given this distribution, multi-particle 
production is determined event-by-event from the CYM 
equations. 

By way of contrast, for instance in the MC-Glauber 
model of Ref. [2f|[26[, a Gaussian distributed energy den- 
sity is added for each participant nuclcon. Its parame- 
ters are the same for every nucleon in every event, with 
the width tuned to be 0.4 fm to best describe anisotropic 
flow data. Unlike the IP-Glasma model, the MC-KLN 
models employ k±_ factorization approximations to de- 
scribe gluon production on the average. Further, as noted 
in Ref. (2?j|, they do not include NBD fluctuations on 
the scale 1/Q S - Finally, in distinction to all the stated 
models, the IP-Glasma model does not assume instan- 
taneous thcrmalization; it includes the effects of pre- 
cquilibrium flow through the Yang-Mills evolution until 
the time where the components of the stress-energy ten- 
sor are matched to those of viscous hydrodynamics. 

Though the IP-Glasma approach is a significant im- 
provement relative to other models with regard to includ- 
ing pre-equilibrium flow, we emphasize that it does not 
include all essential features of the pre-equilibrium stage. 
These include quantum fluctuation driven instabilities 
that change the character of the CYM flow very signifi- 
cantly already by times r ~ l/Q s , and can generate sig- 
nificant amounts of flow. Work in the direction of includ- 
ing these initial quantum fluctuations is highly advanced 
already [28| and "proof of concept" studies demonstrat- 
ing their role in isotropization of 4 theories with heavy- 
ion like initial conditions have been performed |29l - t3l| . 
It is out of the scope of this work to address these pre- 
equilibrium flow studies but we anticipate including them 
in future work. 

The paper is organized as follows. In Section 2, we 
discuss the IP-Glasma model which combines the IP-Sat 
model of nuclcon gluon distributions with classical Yang- 
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Mills evolution of a nuclear collision. Color charge distri- 
butions are sampled from a Gaussian distribution with 
variance g 2 p\(x±). This variance is the color charge 
squared per unit area constructed from the individual nu- 
cleon g 2 /! 2 and is proportional to Qlg A (xj_). Given the 
color charge distribution, the classical Yang-Mills fields 
for each nucleus can be determined, and a unique solu- 
tion for the gauge fields immediately after each collision 
obtained. These then provide the initial conditions for 
numerical solutions of the Yang-Mills equations. The 
details of the numerical computation are discussed in 
Section 3. In Section 4, we discuss results for the fol- 
lowing: i) the dependence of the initial multiplicity as a 
function of the number of participants iVp ar t compared 
to data at RHIC and LHC energies, ii) the multiplicity 
distributions, and the comparison of these to negative 
binomial distributions whose parameters have the form 
conjectured in the Glasma flux tube picture, iii) even and 
odd eccentricity moments and their comparison to those 
in the MC-KLN model. In the final section, we out- 
line the continuation of this work (in a follow up paper) 
where event- by- event hydrodynamic equations are solved 
with input from the IP-Glasma model to determine the 
moments of flow distributions. The IP-Sat model is dis- 
cussed in an appendix. 



II. THE IP-GLASMA MODEL 

We will present in this section details of the IP-Glasma 
model first introduced in Q. The present discussion ad- 
ditionally includes time evolution with the CYM equa- 
tions in 2+1 dimensions. The IP-Glasma framework uses 
the IP-Sat model 0, 0| to determine fluctuating configu- 
rations of color charges in two incoming highly energetic 
nuclei. A review of the IP-Sat model can be found in 
Appendix [5] 

The first step in this framework is to determine the 
event- by- event nuclear color charge densities of each of 
the nuclei which generate the corresponding classical 
fields. This is done as follows. Nucleon positions in the 
transverse plane of both nuclei arc sampled from a Fermi 
distribution 



«(r) = k q 



){r/Rf 



l + exp(^) 



(3) 



where kq is the nucleon density, R the nuclear radius, a 
the skin depth and w denotes deviations from a spherical 
shap<||. 

We next determine the color charge squared per unit 
area of both nuclei g 2 p, 2 A ^ B ^{x , xj_) by summing the cor- 
responding quantities g 2 p 2 (x, bj_, xj_) of all individual 



For 197 Au, the parameters are, R = 6.38 fm, a = 0.535 fm, 
w = 0, and for 207 Pb, R = 6.62 fm, a = 0.546 fm, w = 0. The 
overall normalization kq is irrelevant for sampling the positions. 



nucleons in the nucleus. Here x = (p±)/y/s for zero ra- 
pidity, where (p±) is the average transverse momentum 
of charged hadrons in p+p collisions at a given y/s and 
xj_ is the coordinate in the plane transverse to the beam 
line. For g 2 p, 2 (x, b^, x^) it determines the position of 
the nucleon's center, while is the impact parameter 
relative to its center. The g 2 p, 2 (x,h±) are in turn pro- 
portional to the saturation scale Q 2 ^(x,h±) provided 
by the IP-Sat dipole cross section for each nucleon as 
discussed in Appendix |Al The exact numerical factor 
between g 2 p?(x,h±) and Q 2 ( p - ) (x,h±) depends on the 
details of the calculation (32|; we will discuss it further 
below. 

The distributions g 2 n 2 A ^ B ^(x, xj_) are shown in Fig.[T] 
The degree of lumpiness in this quantity is determined by 
the nucleon size. Given g 2 n 2 A ^(x, xj_), one can sample 




FIG. 1. The incoming color charge densities g^A(B) for two 
gold nuclei at 
to red. 



200 GeV. Increasing density from yellow 



/5 a (xj_) in each event from the Gaussian distribution 

(pV)(xx)pWy±)) = g 2 » 2 A { B)(^±)$ ab s (2 Hx±-y±) ■ 

(4) 

The Gaussian sampled color charges p a (x T , x j_ ) in the 
IP-Sat/MV model act as local sources for small x classi- 
cal gluon Glasma fields. Here, a is a color index running 
from 1 to N 2 — 1. We introduced a general dependence 
on x~ or x + , which depending on the direction the nu- 
cleus is moving, is the longitudinal spatial light-cone co- 
ordinate. Generally, light-cone quantities arc defined as 
v ± = (v° ± w 3 )/-\/2, which translates to proper time and 
spatial rapidity as r = V2x + x~ and r\ = 0.51n(x + /a;~). 
This dependence allows for a finite longitudinal width of 
the nucleus, which we will use in the numerical calcula- 
tion below. 

The classical gluon fields are determined by solving the 
classical Yang-Mills equations 



with the color current 



(5) 



(G) 
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generated by a nucleus A (B) moving along the x + (x~) 
direction (the upper index is for nucleus A). In © we 
have assumed that we are in a gauge where A^ = 
such that temporal Wilson lines along the x + (x~) axis 
become trivial unit matrices. 

The solution to Eq. ([5]) is most easily found in Lorentz 
gauge d^A^ — 0, where the equation becomes a two- 
dimensional Poisson equation 



whose solution can formally be written as 



-PA(B)(x T ,X±)/V 2 L . 



(7) 



(8) 



It will be more convenient to work in light-cone gauge 
A + (A~) = when computing the gluon holds after the 
collision. The solution in this gauge is obtained by gauge 
transforming the result in Lorentz gauge using the path- 
ordered exponential 



A(B) 



'xjj = Pcxp 



-IfJ 



dx 



^A(B) 



(■<- 



(9) 



giving the pure gauge holds [10|, [33J, |3 



A^ (B) (x x ) = e{x-(x + ))-v A[B) {^)d^l [B) {^) , (io) 



A~(A + ) = 0. 



(11) 



The infrared regulator m in Eq. (0) is of order Aqcd and 
incorporates color confinement at the nucleon level. 
Physically, the solution (jlOllip is a gauge transform of 
the vacuum on one side of the light-cone and another 
gauge transform of the vacuum on the other side. We 
have chosen one of them to be zero as an overall gauge 
choice. The discontinuity in the fields on the light-cone 
corresponds to the localized valence charge source [f|. 

The initial condition for a heavy-ion collision at time 
t = is given by the solution of the CYM equations 
in Fock-Schwinger gauge A T = (x + A~ + x~ A + )/t = 0, 
which is a natural choice because it interpolates between 
the light-cone gauge conditions of the incoming nuclei. It 
is also necessary for the Hamiltonian formulation that we 
adopt (gauge links in the temporal (r) direction become 
unit matrices in this gauge). It has a simple expression 



4 Other prescriptions which do not explicitly introduce a mass |35| 
are feasible but they all inevitably involve introducing a nucleon 
size scale. This is because there is a Coulomb problem in QCD 
which is cured only by confinement. The presumption here is 
that physics at high energies is dominated by momenta ~ Q s 
and is insensitive to infrared physics at the scale m. From a 
practical point of view, we observe that our results are weakly 
sensitive to small variations in the scale m. 
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in terms of the gauge fields of the colliding nuclei rl|5l.l36|: 

(12) 

(13) 

(14) 
(15) 

In the limit r -> 0, A n = -E n /2, with E v the longitu- 
dinal component of the electric field. At r = 0, the only 
non-zero components of the field strength tensor are the 
longitudinal magnetic and electric fields, which can be 
computed non-perturbativcly. They determine the en- 
ergy density of the Glasma at r — at each transverse 
position in a single event @, . 

The Glasma fields are then evolved in time r accord- 
ing to Eq. ([5]). Over a time scale ~ 1/Q S the fields are 
strong and the system is strongly interacting. Due to the 
expansion of the system, the fields become weak after 
this time scale and the system begins to stream freely. 
Incorporation of quantum fluctuations in a 3+1 dimen- 
sional CYM simulation will however lead to instabilities, 
which will modify this behavior and potentially keep the 
system strongly interacting for a more extended period 
of time [13, UK. As noted previously, these instabilities 
could isotropize the system, naturally leading to a tran- 
sition to viscous hydrodynamic behavior. The detailed 
study of instabilities and the origin of isotropization is 
a complex task and beyond the scope of this work. For 
recent progress in this direction see [28l . I39rj4l| . We em- 
phasize that key aspects of this work, the event- by- event 
determination of color charge distributions and solutions 
of Yang-Mills equations will be essential ingredients in 
these generalized frameworks as well. In particular, in 
the framework of Ref. pcf . the additional ingredient is 
repeated solution of the CYM equations with slightly 
different seeds drawn from an initial spectrum of fluc- 
tuations. 

III. NUMERICAL COMPUTATION 

We will now discuss the numerical implementation of 
the continuum discussion in the previous section. Be- 
cause the classical gauge field configurations are boost 
invariant, our computations are carried out on 2+1- 
dimensional lattices. From the nuclear color charge den- 
sity squared, determined as described in the previous sec- 
tion, we can sample independent color charges p a (xj_) 
(suppressing x from now on) according to 



.9 2 /4( x -0 



(16) 



The metric in the (r, Xj_,r]) coordinate 
diag(l, -1, -1, -t 2 ) so that A v = -r 2 A' ) . 
of the gauge field are related by 



±x ± A r >. 



system is g M „ = 
The ± components 
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where the indices k, I 
cretized x~ coordinate 



l,2,...,iV y represent a dis- 
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. We typically use N y 
and to = Aqcd ~ 0.2 GeV and will use a value of 
Q s ~ 0.75 g 2 [i, which is somewhat larger than the value 
determined in [32| for the adjoint saturation scale. There 
is some uncertainty in this relation and we have chosen 
0.75 for comparison to experimental data. Changing this 
value primarily changes the overall normalization of the 
multiplicity and energy density. 

For large nuclei, the use of local Gaussian color charge 
distributions is a valid approximation [ltj HH, [42j ■ Mod- 
ifications to Gaussian distributions, relevant for smaller 
nuclei, have recently been explored in (43|. 

On the lattice we work with link variables U^ B ^ - in- 
stead of explicit gauge fields. To determine ^Vb) j we 
first need to compute the path ordered Wilson line ^ in 
its discrctized version 



Va(b)(x±) = Y[ exp 



A(B), N~ 

Pk ( X -U 



(17) 



fc=i 



This quantity already allows us to demonstrate the de- 
gree of correlation and fluctuations in the gluon fields 
of the incoming nuclei. Fig. [5] shows the correlator 
Re[Tr(v], (B) (0,0)V A(B) (x,y))]/N c for = 2760 GeV. 
The characteristic correlation length is 1/Q S , leading to a 
finer granularity El than the nucleon size scale (cf. Fig.[TJ. 
See also (4|| where the evolution of this structure with x 
was computed by solving the JIMWLK renormalization 
group equations. 




FIG. 2. The correlator l/iV c Re[Tr(V^ CB) (0, 0)V A (.b)(x, y))) 
showing the degree of correlations in the gluon fields for lead 
ions at Vs = 2760 GeV. 

To each lattice site j we then assign two SU(iV c ) matri- 
ces V(A),j and V(B),j, each of which defines a pure gauge 
configuration with the link variables 



U A(B),j ~ V A(B),jV A{B) 



(18) 



where +ej indicates a shift from j by one lattice site 
in the i = 1,2 transverse direction. Note that Eq. (TTg|) 



This energy dep endence of the granularity is opposite to the one 
modeled in [Zil . 



is a gauge transform of the unit matrix, hence a gauge 
transform of the vacuum A 1 = 0. 

The link variables in the future light-cone [/*, are de- 
termined from solutions of the lattice CYM equations at 
r = 0, 



tr 



U, 



(A) 



u, 



(B 



(1 



(l + tf*)(c/* +£$,)]}=<), 



U ir ) 



(19) 



where t a are the generators of SU (N c ) in the fundamental 
representation (The cell index j is omitted here) . Eq. (fT9|) 
reduces to Eq. (|12j) in the continuum limit, which can be 
shown by expanding all links for small a: U l ~ l + iagA*. 
The N% — 1 equations (fT9|) are highly non-linear and for 
N c = 3 arc solved iteratively. 

We further need the lattice expression corresponding 
to Eq. (fl3|) . which is the longitudinal electric field in the 
forward light-cone Q 



T Q E [ ( U (A)(^) C/ ( J s) (xx)) (t^t( Xx ) - 1) 
h.c- (^(xx-ir)-!) 



i=x,y 



(13)1 



(20) 



where we have indicated the cell by its coordinate xj_ 
instead of j for clarity. — iy indicates the shift in the —i 
direction by one lattice spacing. 

The lattice Hamiltonian in the boost invariant case is 
given by 



— tr f.'/.' + ^-(N c - Retr U x 2 ) 
t g z a 



-tr 



- 2 + t£ 



tr 



(21) 



where the sum is over all cells in the transverse plane. 
For clarity, we have omitted the cell index j for all quan- 
tities in this expression. <fi is a scalar field, resulting from 
when disallowing rj dependent gauge transformations, 
and 7T = E v = 4>/t is the longitudinal electric field. E l 
with i G {1,2} are the components of the transverse elec- 
tric field that initially are zero. The parallel transported 
scalar field in cell j is given by 



and the plaquette is given by 
ri jf2 



Ui 



] - U 1 !/ 2 - f/ lf f/ 2t 



(22) 



(23) 



With only longitudinal fields present after the collision, 
the total energy density on the lattice at r = is given 
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by 



IV. RESULTS 



e(r = 0) = 



9 2 a 4 



(JV c -Retr0i >2 ) 



(24) 



where the first term is the longitudinal magnetic, the 
second the longitudinal electric energy density. 

To evolve the system from t = forward, we use the 
Hamiltonian equations of motion obtained from taking 
the Poisson brackets of the fields with the Hamiltonian 
(ED 



Ui = i—E l Ui (no sum over i) 

T 

(f> = TTT 

E x = ^Pi,2 + Ui,- 2 - U{ 2 - U{_ 2 - Ti] 



+ -[01,0] 
T 

■ ^ [C/2,1 + C/2,-1 - ^,1 - 4-1 - T 2 ] 

+ -[01,0] 
T 



E 2 



+ 4>-i - 20 



(25) 
(26) 

(27) 

(28) 
(29) 



where T x = ^tr[U h2 + U\,-i - U\ >2 - U\_ 2 ) 1, and 

T 2 = ^tr[[C/ 2jl + U 2 ,-i - UI A - Ul_ x ] 1 with the N c x 
N c unit matrix 1. The subtraction of these traces takes 
care of the fact that the components of F^ v are traceless 
color matrices. Note that the Hermitian conjugates of a 
plaquette simply reverse direction. 

To determine the energy density at times r > we 
evaluate the rr-component of the energy-momentum ten- 
sor at site j 



rpTT 

3 



J^tr [ E l,j + E lj+e 2 + E y,j + E l,j+eA 



2r 2 



rtr 



+ {4>3-e^-U x , j -^jUl i _ gi f 



2r 2 



tr 



+ (<p3-i 2 -Uy, j -^ j Ul j _ g2 f 
4D y 



(30) 



where YIad indicates the sum over the plaquette defined 
in (|23|) and the three other plaquettes with the same 
orientation connected to site j. Note that care has to be 
taken that all quantities are defined at the same point, 
in this case lattice site j. 



A. Energy density and dN g /dy 

Having outlined the numerical procedure, we will now 
present results from the CYM evolution of the IP-Sat 
initial conditions. We begin by showing the evolution of 
the energy density (|30|) . integrated over the transverse 
plane in Fig. [3] As can be seen from Eq. (|2"4")l , ini- 
tially the system consists of only longitudinal chromo- 
magnetic and chromo-electric fields. Over the course of 
the evolution, the transverse components of the energy 
density grow and catch up with the longitudinal modes. 
The expansion of the system leads dE/rdy to drop like 
1/t and dE/dy becomes a constant. 
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FIG. 



3. Components of the integrated energy density dE/rdy 
as a function of r from one single event. Labels indicate 
components of dE/rdy. 

Fig.S] shows the structure of the energy density in the 
transverse plane after the collision, at time r = Ofm 
(Eq. (|24[0 and after classical Yang- Mills evolution in 2+1 
dimensions for At = 0.2 fm, which is of the order of 
1/Qs, the time scale over which interactions are still sig- 
nificant. After this time, expansion causes the fields to 
become weak and the system becomes freely streaming. 
As shown in Fig. 02 at the later time all four components 
of Eq. (|2"Tj) contribute. As can be expected, the distribu- 
tion becomes smoother with the evolution. 

We next compute the gluon multiplicity per unit rapid- 
ity dNg/dy. Since the multiplicity is not a gauge invari- 
ant quantity, we fix transverse Coulomb gauge {diA 1 = 0, 
with i summed over 1, 2). Then the lattice expression for 
dN g /dy is given by [|,|4| 



dNg 

dy 



+ rtr(7r(k x )7r(-k ± )) , (31) 



assuming a free massless lattice dispersion relation for 
the interacting theory. This leads to the appearance of 
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FIG. 4. (Color online) Energy density (arbitrary units) in the 
transverse plane at r = Ofm (upper panel) and r = 0.2 fm 
(lower panel). The structures are smoothed by the evolution 
over the first At ~ 1/Q S . 



the square root of 



2 k y 



(32) 



The result for (dN g /dy)/(N pllIt /2) vs. iV part at time 
t = 0.4 fm. where iV part is the number of participant 
nucleons, is shown in Fig. [5] The IP-Glasma model does 
not have the concept of "wounded nucleons" because we 
treat the nucleus as a coherent system of gluon fields 
correlated on distance scales 1/Q S much smaller than the 
size of a nucleon. To determine iV par t we use the same 
Monte Carlo Glauber (MC-Glauber) model as employed 
by the experimental collaborations. Nucleons that were 
sampled for each nucleus, as described in Section IIIIl 
are assumed to be participant nucleons if the relative 
transverse distance between them and a nucleon from the 
other nucleus is smaller than D = \J a nn/^, where <jnn 
is the total inelastic cross section, (Jnn = 42 mb for y/s = 
200 GeV Au+Au collisions at RHIC and 64 mb for yjl = 
2.76 TeV Pb+Pb collisions at LHC. The reader should 
note that the MC-Glauber model and the inelastic cross 
sections are solely used to determine V par t to compare 
to experimental data. They are not an input for the IP- 
Glasma model. 

In the upper panel of Fig.O we show the multiplic- 
ity distribution as a function of iV part for fixed coupling. 
We multiplied the gluon multiplicity by 2/3 to convert to 
charged particle multiplicity. Note that the overall nor- 
malization is chosen (by varying the ratio between Q s 
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FIG. 5. Gluon multiplicity (dN g /dy)/(N ps . It /2) at r = 
0.4fm/c times 2/3 compared to experimental charged par- 
ticle (dN/dy)/(N pa _vt/2) data for yjl = 200 GeV Au+Au and 
\fs — 2.76 TeV Pb+Pb collisions as a function of iV pa rt for 
fixed coupling (upper panel) and running coupling (lower 
panel). The pale blue and red bands are a collection of the 
multiplicities for individual events, with the solid lines repre- 
senting the average multiplicity. Experimental data from [47l ] 
and Eil. 



and g 2 /i and a s in the fixed coupling case) to agree with 
the RHIC data for charged particles. This allows us to 
better compare the shape of the result and the experi- 
mental data. (The pale bands denote results from the 
individual events and demonstrate the range of fluctu- 
ations around the mean. See below for more details.) 
However, we know that there is entropy production in 
the system and the initial gluon multiplicity should not 
account for all observed final particles. The logarithmic 
uncertainty in Q s as well as some numerical uncertainty 
(for details see [THHU) in the factor between Q s and g 2 [i 
introduce some freedom that allows to adjust the normal- 
ization of the initial dN g /dy. This also allows to adjust 
the energy density when fine tuning to experimental data 
when using this model with a viscous hydrodynamic evo- 
lution model that accounts for entropy production. 

While the RHIC result is reasonably well described, 
both the normalization and shape of the LHC result dis- 
agree strongly with the experimental data for the fixed 



8 



coupling case. Note that the IP-Sat model naturally gen- 
erates different x dependencies of Q s for protons and nu- 
clei, the latter giving a stronger dependence on x. Be- 
cause we use the parametrization of Q s from DIS off of 
protons to construct nuclei, we get the right spatial de- 
pendence of color charge distributions but the x depen- 
dence for nuclei is off, and the stronger a;-dcpcndcnce has 
to be introduced by hand. We observe that our results 
are consistent with Ref. [l8j . 

The inclusion of running coupling as shown in the 
lower panel of Fig. [5] improves the agreement with the 
LHC data significantly. The shape is now very well de- 
scribed apart from a very modest under-prediction at 
small iVpart. This discrepancy could be from entropy 
production in the later evolution stages, whose relative 
contribution can be expected to be somewhat larger for 
smaller systems. The scale ft for the running of 



a a {ft) 



2tt 



(ll-fJV»la(/i/AQCD) 



(33) 



was chosen to be ft = (max(Qf(x±),Qf(x±))), where 
the average is over the whole transverse plane. As usual, 
there is some ambiguity in choosing the scale, but we have 
verified that using the minimum instead of the maximum 
or twice the value did not make a significant difference 
for the shape of the result. Note that the CYM evo- 
lution does not depend on a s as it scales out from the 
equations of motion and only enters in the final compu- 
tation of the multiplicity. Furthermore note that the use 
of MC-Glauber to determine N palt is crucial for repro- 
ducing the experimentally observed shape of the result. 
In this way, we include both fluctuations in iV par t (along 
the horizontal axis) as well as in dN/dy (along the ver- 
tical axis). Individual events (15000 each) are plotted as 
scattered points in the background to show the degree of 
fluctuations explicitly. In a forthcoming paper, we will 
evaluate how event- by- event viscous hydrodynamic sim- 
ulations modify this computation. 



B. Multiplicity distributions 

In Fig. [5] we present the probability distribution of 
dNg/dy at RHIC energies. An essential ingredient is the 
probability distribution of impact parameters, which is 
determined by the Glauber model to be 



dP, 



incl 



1 - (1 - <?nnTab) 



AD 



d 2 b 



j d*h x {l- {I - <t nn T ab )AB)> 



(34) 



with the overlap function Tab- One could in principle 
compute this distribution in the Glasma framework, but 
a first principles computation is extremely difficult. The 
probability for no inelastic interaction (an essential in- 
gredient in the above equation) at a given impact pa- 
rameter requires an understanding of diffractive/elastic 
interactions which is incomplete at present in all QCD 
based frameworks. The Glauber model, with parameters 



tuned to data, is therefore a good effective model for this 
aspect of our computation. 

We compute the n-particle multiplicity distribution by 
first sampling the impact parameter b from a uniform 
distribution, computing the resulting dN g /dy from the 
IP-Glasma model, and when binning these values into 
the histogram shown in Fig.[6j weight the result with a 
factor 27r b dP\ nc \/d bj_ depending on the b value used in a 
given event. The STAR data shown is uncorrected, which 
makes a direct comparison difficult. We therefore scale 
dN/dy by a factor of 0.8 and achieve good agreement 
with the data. Note, however, that the scaling factor 
should depend on dN/dy (49j . 



RHIC 200 GeV 
STAR (uncorrected) — IP-Glasma 

— IP-Glasma (dN/dy x 0.8) I 




400 600 

dN/dy 

FIG. 6. Probability distribution of gluon multiplicities dN/dy 
at r = 0.4fm/c. Shown are also the distributions for some 
limited ranges of impact parameter 6, which are described 
by negative binomial distributions. Experimental data from 
STAR H. 




1000 1500 

dN/dy 

FIG. 7. Same as Fig. ©for LHC energy y/a = 2.76 TeV. 

We also show three distributions obtained by con- 
straining the impact parameter range to demonstrate 
that their shape resembles a negative binomial distribu- 
tion (NBD). The negative binomial fluctuations of the 
transverse energy as emerging from the IP-Glasma model 
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have previously been discussed in Q . Here we show that 
they also appear in the multi plic ity distribution. In the 
Glasma flux tube framework [50 . l51| , negative binomial 
distributions 



T(k + n) n n k k 
T(k)T(n + l) (n + k) n+k ' 



with 



(35) 



(36) 



arise |15| . with k inversely proportional to the width of 
the NBD. k here is proportional to the number of flux 
tubes QgSj_, where S± is the transverse size of the sys- 
tem. (Hence smaller systems at a fixed energy generate 
more fluctuations.) In the expression above, £ is an in- 
trinsically non-pcrturbativc function which can be com- 
puted ab initio from solutions of the CYM equations. To 
determine £ it is sufficient to compute the double inclu- 
sive distribution by solving CYM equations as done in 
Ref. fl6j . Here we are computing the n-particle inclu- 
sive distribution, and can thus also extract £ if the NBD 
description of these distributions is robust. £ in general 
can depend on Q 2 S±; however, a powerful test of how 
robust the Glasma flux tube picture is depends on this 
dependence being weak for large values of Q 2 S±. Wc 
will therefore discuss the behavior of £ below and what 
the results tell us about the nature of different sources of 
quantum fluctuations. 

Within the IP-Glasma model, we determine ((Q 2 S±) 
by extracting k from a fit with an NBD and computing 
an average (Q 2 S±) by summing over the minimum of the 
two Q 2 in the whole transverse plane. 

Wc first determine n and k of Eq. (|35p from the fit to 
the multiplicity distributions at fixed impact parameter 
b. The results are shown in Fig. [8] Interestingly, the 
ratio k/n is greater than one for central collisions and 
at large impact parameters approaches approximately 
k/n 0.14, which is close to the value determined from 
fits to distributions in p+p collisions [27| . 

The corresponding £- values are shown in FigJH] as a 
function of the average values of Q 2 S± for a given b. Wc 
observe a strong dependence of £ on Q 2 S±, which is in 
disagreement with the flux tube picture. The reason is 
that the effect of fluctuations in the number of wounded 
nucleons (which were not considered in the derivation of 
Eq. (|36|) ) in addition to fluctuations in the color charge 
distributions, make the distribution wider (£ smaller), es- 
pecially at large impact parameters (small Q 2 S±) where 
geometrical fluctuations dominate. 

This behavior of £ is compatible with previously ex- 
tracted values from fits based on k± factorization to 
A+A multiplicity distributions [HI, where small dN/dy 
required small values of zeta and large dN / dy larger val- 
ues to achieve a good fit. The IP-Glasma model auto- 
matically produces this variation of £, leading to very 
good agreement with the experimental data as shown in 
Fig.® 
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FIG. 8. Average gluon number n and width parameter k as 
a function of impact parameter b at ^/s = 200 GeV. 
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FIG. 9. Proportionality factor £ at RHIC energy in k — 
(£(iV c 2 - l)/(2n))Q 2 a S± as a function of Q*S± including all 
fluctuations. QlS±_ corresponds to the average value for a 
given b. b was varied in 1 fm steps from to 13 fm. 



To be able to better compare to the flux tube picture 
and Eq. ([551 , we now consider only the effect of fluctu- 
ations in color charges. To do so we average over the 
nuclear color charge density squared over many nucleon 
configurations. This results in a smooth distribution that 
removes fluctuations in the wounded nucleon number and 
positions. 

This is in the spirit of the original Glasma flux tube 
perturbative [l5T | and non-perturbative [l6[ computa- 
tions. The result for this £ S mooth is shown in Fig. 1101 
After starting out near 1, £ sm ooth drops and approaches 
a constant value of approximately £ sm ooth = 0.2 for large 
Q 2 S±. This means that at low parton densities, k is ini- 
tially more constant than expected in the Glasma flux 
tube picture but then becomes proportional to Q 2 S± at 
high parton densities as anticipated. 

The fact that £ at small impact parameters (cf. Fig.|H]) 
approaches that in the smooth case for the same Q 2 S± 
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FIG. 10. Proportionality factor ( in k = (CC-ty? - 
1) / (2ty))Q 2 S± as a function of Q 2 S± for averaged nucleon po- 
sitions (squares) and with nucleon fluctuations at fixed impact 
parameter b — Ofm (circles). At large Q 2 S± the result for the 
smooth distribution approaches a constant as predicted by 
the Glasma flux tube model for n-gluon correlations. The re- 
sult for fluctuating nucleon positions at constant b = Ofm is 
very similar and becomes very weakly dependent on Q 2 S±. 



gives hope that there is a chance to experimentally ac- 
cess a regime where the flux tube picture is valid. Fixing 
b = fm and increasing the energy, and this way in- 
creasing Q 2 S S± while keeping S± as constant as possible, 
reduces fluctuations in the nucleon number. Indeed we 
find that the result for the extracted (, shown as blue 
circles in Fig.[T0l is very close to the one obtained with 
smooth initial distributions, and its dependence on Q^S± 
becomes weak at large QlS±. 



C. Eccentricities 

In Ref . Q , we presented results for £2 and £3 defined in 
(TJJ and compared to results from an MC-Glauber model 
and the MC-KLN model [H, [H [H|. Here we extend 
this study by comparing eccentricities up to e§. The re- 
sults are shown in Fig. [IT] and the conclusions that can 
be drawn are mainly that the purely fluctuation driven 
odd harmonics £3 and £5 from the IP-Glasma model are 
larger than those from the MC-KLN model [54[ for all b, 
while £2 is smaller than that computed in the MC-KLN 
model, in particular for b > 3 fm. As a consequence, the 
ratio £2 1 £3 is smaller than in the MC-KLN model, which 
is going to decrease the ratio of V2 /«3 obtained after hy- 
drodynamic evolution, making it more compatible with 
experimental observation. £4 and £§ are almost equal or 
larger than those from the MC-KLN model. We make the 
comparison to MC-KLN at t = Ofm/c, because at later 
times we would also have to take into account the pre- 
equilibrium flow built up in the CYM simulation. This 
effect will be included in detailed event- by- event simula- 



tions that convert the spatial anisotropics into momen- 
tum anisotropics in a follow up paper in this series. To 
show the effect of the CYM evolution on the eccentricities 
themselves we show as an example their time evolution 
for b = 8fm in Fig. [12] As expected already in 0], the 
change in all e n is very weak over the first 0.4 fm/c. After 
this time all £ n begin dropping as the systems is freely 
streaming and hence becoming more isotropic. 
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FIG. 11. Even (upper panel) and odd (lower panel) eccen- 
tricities from the IP-Glasma model compared to those from 
MC-KLN. 



V. SUMMARY AND OUTLOOK 

This paper expands on the IP-Glasma model of fluctu- 
ating initial conditions for heavy-ion collisions first pre- 
sented in Ref. 0- More details of the computations are 
presented as well as novel results for the single inclu- 
sive and n-gluon multiplicity distributions and higher 
even and odd eccentricity moments. In addition, the 
CYM equations are solved for finite proper times r unlike 
Ref. Q where results were extracted only for r = 0. 

We observed that the running coupling CYM results 
give good agreement with the RHIC and LHC single in- 
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FIG. 12. Time evolution of the eccentricities from the IP- 
Glasma model at impact parameter b = 8fm. 



elusive multiplicity distributions as a function of N pax t. 
The computed n-particle inclusive multiplicity distribu- 
tion shows good agreement with the uncorrected STAR 
data on the same once a constant correction factor is ap- 
plied. A prediction is made for multiplicity distributions 
at the LHC. We further observe that our results are well 
described as a convolution of negative binomial distribu- 
tions at different impact parameters. The parameters of 
the negative binomial distributions are extracted, and it 
is observed that, in the approximation of smoothed nu- 
cleon configurations, for large parton densities, the pre- 
dictions of the Glasma flux tube picture are recovered. 
For the realistic situation of fluctuating wounded nucleon 
configurations, one still obtains NBDs albeit with signif- 
icantly wider widths. The non-perturbative coefficient £ 
introduced in the Glasma flux tube description quanti- 
fies the effect of wounded nucleon fluctuations. Predic- 
tive power, in particular for central impact parameters 
at higher energies, is still retained because £ is (nearly) 
energy independent, while the width parameter k (in 
Eq. [36]) has a strong energy dependence controlled by 
the saturation scale Q 2 . 

The computation of eccentricity moments is performed 
up to eg. It is seen that £2 is smaller than in the MC- 
KLN model (used as an initial condition in many hy- 
drodynamic studies), while the odd moments are larger, 
pointing to the additional role of multiplicity fluctuations 
in the IP-Glasma model. 

An essential follow up to this work is to match the re- 
sults of the IP-Glasma model, event- by- event, to viscous 
hydrodynamic simulations. This will allow one to gauge 
the effects of dissipative flow in modifying the energy 
and multiplicity distributions and on the conversion of 
spatial anisotropics into momentum anisotropics. These 
have the potential to significantly enhance our under- 
standing of the transport properties of the quark-gluon 
plasma, with the caveat that a systematic treatment of 
pre-equilibrium flow including instabilities can alter some 
of these conclusions significantly. 
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Appendix A: The IP-Sat model 

The impact parameter dependent dipole saturation 
model (IP-Sat) [J] is a refinement of the Golec-Biernat- 
Wiisthoff dipole model [5^, [5(1 to give the right pertur- 
bative limit when rx — >• 3]. It is equivalent to the 
expression derived in the classical effective theory of the 
CGC, to leading logarithmic accuracy (57j . 

The proton dipole cross-section in this model is ex- 
pressed as 



dcr gj P 
d 2 b^ 



(rj., x, b x ) = 2 [1 - exp (-F(r x ,x, bx))] (Al) 



with 



F(rx,z,bx) = ^r ± 'a s ((i')xg(x,tf)T p (bx).(A2) 



Here the scale p? is related to dipole radius rx as 



A 2 



- 2 
Mo 



(A3) 



and the leading order expression for the running cou- 
pling is given by Eq. (|33|) . The model includes saturation 
as eikonalized power corrections to the DGLAP leading 
twist expression and may be valid in the regime where 
logs in Q 2 dominate logs in x. The saturation scale for 
a fixed impact parameter is determined self-consistently 
by requiring that the dipole amplitude (within brackets 
in eq. lAl[) have the magnitude N(x, rs, bx) = 1 — e -1 / 2 , 
with Q 2 p = 2/r|. We note that there is an overall loga- 
rithmic uncertainty in the determination of Q 2 p (x, bx)- 
For each value of the dipole radius, the gluon density 
xg(x,jj, 2 ) is evolved from fi 2 , to jx 2 using LO DGLAP 
evolution equation without quarks, 



dxg{x,fi 2 ) _ asijj 2 ) 
d log jl 2 2ir 



1 

J dzP gg {z)- z g[- zl fi 2 ) (A4) 



Here the gluon splitting function with Nf flavors and 
CU=3 and Tr=1 is 



-7% — —rr + ~ — - + z (! _ z ) 

(1 - z)+ z 



11 N 



5(1 -z). 



(A5) 
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The initial gluon density at the scale is taken to be of 
the form 

xg(x,ft)=A g x- x »(l-xf*. (A6) 

An important feature of the IP-Sat model is the In- 
dependence of the dipole cross-section, which is intro- 
duced through a gluon density profile function T(b). This 
profile function is normalized to unity and is chosen to 
have the Gaussian form 

r 'CJ-5s;-'(w)' (A7) 



where Bq is a parameter fit to the HERA diffractive data. 
This corresponds to (b 2 ) = the average squared 

gluonic radius of the proton, which corresponds to a very 
small transverse radius ~ 0.5 fm and a three-dimensional 
radius of - 0.61 fm [H|. 
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